An exact representation of the fermion dynamics in terms of Poisson processes 
and its connection with Monte Carlo algorithms 
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We present a simple derivation of a Feynman-Kac type formula to study fermionic systems. In 
this approach the real time or the imaginary time dynamics is expressed in terms of the evolution 
of a collection of Poisson processes. A computer implementation of this formula leads to a family 
of algorithms parametrized by the values of the jump rates of the Poisson processes. From these an 
optimal algorithm can be chosen which coincides with the Green Function Monte Carlo method in 
the limit when the latter becomes exact. 



A crucial issue in quantum Monte Carlo methods is 
the choice of the most convenient stochastic process to 
be used in the simulation of the dynamics of the system. 
This aspect is particularly evident in the case of fermion 
systems due to the anticommutativity of the vari- 

ables involved which for long time have not lent them- 
selves to direct numerical evaluation. In a recent paper 
|^| progress has been made in this direction by provid- 
ing exact probabilistic expressions for quantities involv- 
ing variables belonging to Grassmann or Clifford alge- 
bras. In particular, the real time or the imaginary time 
evolution operator of a Fermi system or a Berezin inte- 
gral can be expressed in terms of an associated stochas- 
tic dynamics of a collection of Poisson processes. This 
approach depends on an older general formula to rep- 
resent probabilistically the solution of a system of or- 
dinary differential equations (ODE) in terms of Poisson 
processes ||. In this paper, we present a simple deriva- 
tion of a similar formula to study fermionic systems, in 
particular, the Hubbard model. However, the fermionic 
nature of the Hamiltonian plays no special role and sim- 
ilar representations can be written down for any system 
described by a finite Hamiltonian matrix. A computer 
implementation of this formula leads to a family of al- 
gorithms parametrized by the values of the jump rates 
of the Poisson processes. For an optimal choice of these 
parameters we obtain an algorithm which coincides with 
the well known Green Function Monte Carlo method in 
the limit when the latter becomes exact 0. In this way 
we provide a theoretical characterization of GFMC. 

Let us consider the Hubbard Hamiltonian 

|A| |A| _^ 
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where A C Z d is a finite d-dimensional lattice with cardi- 



nality |A|, {1,...,|A|} some total ordering of the lattice 
points, and Ci a the usual anticommuting destruction op- 
erators at site i and spin index a. In this paper, we are 
interested in evaluating the matrix elements (n'\e~ Ht \n) 
where n — (nif,n.i|, . . . , «|a|T; "|a i) are the occupation 
numbers taking the values or 1 M] . The total number 
of fermions per spin component is a conserved quantity, 
therefore we consider only configurations n and n' such 

that J2[=i n 'ia = Si=i n i<7 Ior u =U- I n the following 
we shall use the mod 2 addition n(Bn' = (n + n') mod 2. 

Let T = 1 < i < j < |A| : % ^ 0} and |r| its 

cardinality. For simplicity, we will assume that jjy = r/ if 
(i,j) G r and 7,; = 7. By introducing 

Kja{n) = {n@ \ ia © lj<r\4 a c ja + c] a c ia \n) 

= (— 1)"*"^ H n j-1<7 

X [nj a {n ia e 1) - n irT (n ja © 1)] , (2) 
where \ la = (0, . . . , 0, l iff , 0, . . . , 0), and 

|A| 

V(n) = (n\H\n) = 7^n iT n a , (3) 

i=l 

the following representation holds 

(n'\e- Ht \n)=E(S n ,. n(BNt M t ) (4) 

M l = cxp( E / l0 S [W~%<> © K s )] dN! ja 

-J V(n®N s )ds + 2\r\pty (5) 

Here, {Njj a }, € T, is a family of 2|T| indepen- 

dent Poisson processes with parameter p and TV* = 
(NlpN^, . . . , Nh.pNL,^) are 2|A| stochastic processes 
defined as 
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TV*. 
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(i.j)^T: i—k or j—k 



We remind that a Poisson process AT* with parameter p 
is a jump process characterized by the following proba- 
bilities: 



P (N t+S — AT* = k) = ^V ps . 
v ' kl 



(7) 



Its trajectories are piecewise-constant increasing integer- 
valued functions continuous from the left. The stochastic 
integral J dN* is just an ordinary Stieltjes integral 



f(s,N s )dN s 



E 

k: s k <t 



where Sk are random jump times having probability den- 
sity p(s) — pe~ ps . Finally, the symbol E(. . .) is the ex- 
pectation of the stochastic functional within braces. We 
emphasize that a similar representation holds for the real 
time matrix elements (n'\e~ lHt \n) . 

Summarizing, we associate to each r)ij / a link con- 
necting the sites i and j and assign to it a pair of Poisson 
processes with a =11. Then, we assign to each site 
i and spin component a a stochastic process 7V* CT which 
is the sum of all the processes associated with the links 
incoming at that site and having the same spin compo- 
nent. A jump in the link process iVy CT implies a jump 
in both the site processes N* a and Nj a - Equations (|) 
and (H|) are immediately generalizable to non identical 
parameters Tjij and 7j. In this case, it may be convenient 
to use Poisson processes N*j a with different parameters 

Pija • 

We now show that the representation (^j^) follows from 
the general formula to represent probabilistically the so- 
lution of an ODE system || and the expression of the 
matrix elements of H. The matrix elements (n'\e~ Hl \n) 
obey the ODE system 



-HI 



dt 



n)=-J2(n'\H\n"}{n"\e- m \n) 
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One may check that (||||) is indeed solution of (g) by 
applying the rules of stochastic differentiation. We have 

E{5 n , n(BN t+dtM t+dt ) 

= E- E ( IT II <^',™"©<w 

x e -v( n ® N *)dt+2\r\ P dt Sn „ j^m*] (9) 



For the Markov property, the expectation of the factors 
containing the stochastic integrals in the interval [t, t+dt] 
can be taken separately. By expanding each one of them 



over all possible numbers of jumps of the Poisson pro- 
cesses as 



- 5 , „ e° e~ pdt 

+ On',n"eii CT ©lj CT e L 1 e pdt + 

-5 n ,, n „}pdt + 0(dt 2 ) , 
up to order dt we obtain 
E {5 n i >n@ N t+dt M t+dt ) 
S r . 
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>V(n")dt 
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7 r)\ lJry (n")dt 
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Finally, we rewrite this relation as 
AT) 

= E {5 n , n@N t+dtM t+dt ) - E (6 n > ^(BNtM 1 ) 

]{n'\H\n")E (<5„», nffiJ v^') dt. (11) 
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It is clear that the fermionic nature of the Hamilto- 
nian H plays no special role in the above derivation. 
For later use, note that summing ( |ll| ) over n' we have 
dE(M f ) = -E^M^dt, where 

H l = ^(n'|tf|nffiiV'} 

71 ' 

-V E E A u-( n © N*) + V{n@N t ). (12) 



E 

(i,i)era=u 



In order to construct an efficient algorithm for eval- 
uating we start by observing that the functions 
93 N s ) vanish when the occupation numbers 
ni a © N? a and rij a © N? a are equal. We say that for 
a given value of a the link is active at time s if 
Aijo-(n © AT 6 ) ^ 0. We shall see in a moment that only 
active links are relevant. Let us consider how the stochas- 
tic integral in (|^) builds up along a trajectory defined by 
considering the time ordered succession of jumps in the 
family {Nfj a }. The contribution to the stochastic inte- 
gral in the exponent of (^|) at the first jump time of a 
link, for defmiteness suppose that the link iiji with spin 
component a± jumps first at time s±, is 

log [r?p _1 A ilil(7l (n © AT 51 )] 6(t - s x ), 

where N Sl = due the assumed left continuity. There- 
fore, if the link i\j\<7\ was active at time we obtain 
a finite contribution to the stochastic integral otherwise 
we obtain — oo. If s\ > t we have no contribution to the 
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stochastic integral from this trajectory. If s\ < t a second 
jump of a link, suppose i 2 j2 with spin component 02, can 
take place at time s 2 > si and we obtain a contribution 

log [vp- l Kh*An © N")] 6{t - s 2 ). 

The analysis can be repeated by considering an arbitrary 
number of jumps. Of course, when the stochastic integral 
is — 00, which is the case when some A = 0, there is no 
contribution to the expectation. The other integral in (|^) 
is an ordinary integral of a piecewise constant bounded 
function. 

We now describe the algorithm. From the above re- 
marks it is clear that the only trajectories to be consid- 
ered are those associated to the jumps of active links. It 
can be seen that this corresponds to the conservation of 
the total number of fermions per spin component. The 
active links can be determined after each jump by in- 
specting the occupation numbers of the sites in the set 
r according to the rule that the link ij is active for the 
spin component a if rii G + rij a = 1. We start by de- 
termining the active links in the initial configuration n 
assigned at time and make an extraction with uniform 
distribution to decide which of them jumps first. We 
then extract the jump time si according to the proba- 
bility density PaAs) — ^4ipexp(— A\ps) where A\ is the 
number of active links before the first jump takes place. 
The contribution to A4 l at the time of the first jump is 
therefore, up to the factor exp ( — 2|T|pt), 

Tjp-'X^in © AT^)e- l/ ("® iVS1 )^ e -( 2 l r l-^)^ 
x6(t - 8i ) + e -y(^^)t e -(2\r\-A l)p t (si _ ^ 

where exp[— (2|T — Ai)ps] is the probability that the 
2|r| — A\ non active links do not jump in the time inter- 
val s. The contribution of a given trajectory is obtained 
by multiplying the factors corresponding to the different 
jumps until the last jump takes place later than t. For a 
given trajectory we thus have 



k>l 



Xe [A h p-V( n SN^)]( Sk - Sk _ 1 ) _ 
+e [A,p-F(neJV t )](i- Sfc - 1 ) Q( Sk _ t )j _ ( 13) 

Here, At — A{n © N Sk ) is the number of active links 
in the interval (sk-i,Sk\ and so = 0. Note that the ex- 
ponentially increasing factor exp (2\T\pt) in (|) cancels 
out in the final expression of . The analogous expres- 
sion of M for real times is simply obtained by replacing 
t] — > if] and 7 — > ij. 

Let us specialize the algorithm for the calculation of 
the ground state energy Eq. This can be related to the 
matrix elements (Ef) in the following way 



£o = -lim%^^ 



l\„-Ht 
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-Ht 
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(14) 



The denominator in this expression corresponds to eval- 
uating the expectation in without the 5 and is esti- 
mated by 



n e 



-Ht 



n) 



1 K 

K z_^, Jvl p> 



(15) 



where the index p denotes one of the K simulated trajec- 
tories and A4p is the value of Ai f for the pth trajectory. 
The numerator of (14) is estimated similarly to ( |l5| ) with 
Mp replaced by H p M P , where Hp is the value for the pth 



trajectory of Tl 1 given by (|T^). 

The variance of the stochastic process A4* increases 
with t and its distribution is not well estimated if the 
number K of trajectories remains fixed. As an alterna- 
tive to increasing K, one may resort to the reconfigu- 
ration method ||. A simulation with K trajectories is 
performed for a time t but is repeated R times choosing 
randomly the initial configurations among those reached 
in the previous simulation p 11 1. 

In principle, the algorithms parametrized by p are all 
equivalent as (|[^) holds for any choice of the Poisson 
rates. However, since we estimate the expectation values 
with a finite number of trajectories, this may introduce 
a systematic error. In Fig. 1 we show the dependence 
of the error Eq — £;g xact as a function of p in the case 
of a small one-dimensional system which can be exactly 
diagonalized. It is evident that the best performance of 
the algorithm is in correspondence of the natural choice 
p ~ rj independently of the interaction strength. This 
behavior can be understood on the basis of the following 
qualitative argument. The average number of configura- 
tion changes in the time t is Apt where A is the average 
number of active links. This is also the average number of 
terms in the product of ( |lg| ) . In absence of sign problem 
as in the case of Fig. 1, we roughly estimate the r.h.s. of 
(|l5| ) as (i r ip~ 1 ) Apt e < - Ap_y ' t , where V is the average po- 
tential. We see that the derivative with respect to p of 
this expression vanishes for p — 77. It is easily seen that 
this corresponds also to a minimum of the variance of 
M 1 and, therefore, to a minimum of the statistical error. 

The convergence features of the algorithm can be im- 
proved also using the importance sampling method [0 . 
Consider the operator H isospectral to H defined by 
the matrix elements (n'\H\n) = (n'\g) (n'\H \n) (nig) -1 , 
where \g) is a given state. The stochastic representation 
(^-|^) and the corresponding algorithmic implementation 
hold unchanged for the new operator H with the substi- 



tution Ayv(n) 



\ ija (n)(n 8 l«r ® lja\g){n\g)~ 



In 

this way, however, the value of p is not tuned to obtain 
the best performance. For this purpose we have to choose 
p dependent on the indices ija and on the configuration 
n so that Py CT (n) = 77 \(n © l ia © lj a \g){n\g)~ l I. 
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We consider now the connection of our approach with 
the GFMC method. In the GFMC, the ground state of 
the system is filtered out from a given initial state \n) by 
iteratively applying the operator G = 1 — Hil^ 1 . In the 
occupation number representation we can write 



{n'\G\n)=V n >n b r 



where 



\(n'\G\n)\ 



En"l<n"|G|n>| 
is a stochastic matrix and 



(16) 



(17) 



(18) 



a, possibly negative, weight factor. A trajectory is de- 
fined as a series of steps, each one of duration fi , in 
which the configuration changes from n to n' according 
to the stochastic matrix V n 'n- The weight of a trajec- 
tory is Yi{n} bn'n, where n runs over the intermediate 
configurations visited by the trajectory Under our 
hypothesis ?7y = 77 if G T, the GFMC algorithm 

becomes very simple because all the active links have the 
same probability to jump. 

For large values of f2 the jump probability vanishes 
and the GFMC as described above becomes inefficient. 
As remarked in jjj, one can cope with this difficulty in 
the following way. The probability that a configuration 
n changes at the step n s + 1, being unchanged in the pre- 
vious n s steps, is p(n s ) = (V nn ) ne (l - V nn )- Therefore, 
we can directly assign the weight 

{b nn ) n "b n , n = X ija (n) {1 + [ V A(n) - V(n)} fT 1 }™^ 1 

to a configuration change n — > n' = n © 1^ © lj a in 
which a fermion with spin a is moved from the site i 
to the site j or viceversa. As before, A(n) is the num- 
ber of active links in the configuration n and \ija(n) 
and V(n) are given by (||) and (||), respectively. The 
random integer n s is extracted with probability p(n s ), 
e.g. n s = [logr/ logPjmJ , where r is a random num- 
ber with uniform distribution in [0,1]. Note that the 
distribution of n s becomes Poissonian with parameter 
— \ogV nn ~ r]A(n)^ 1 for large. At this level, there 
is no apparent connection between this Poisson process 
and those entering in our formula (^J). However, as sug- 
gested in H one can take the continuum limit SI -1 — > 
and reconstruct in this way the semigroup exp(-Ht). 
It is easy to verify that in this limit the GFMC algo- 
rithm coincides with that obtained from our formula for 
the optimal choice p — 77. In fact, the time interval be- 
fore a configuration change takes place is s = n s f2 -1 . 
For SI -1 — > 0, the random variable s becomes con- 
tinuously distributed in [0, 00) with probability density 
p(s) = Q,p{n s ) = f]A(n) exp[— r]A(n)s] and the trajec- 
tory weight corresponding to the interval s reduces to 



^ij<r(n) exp {[i]A(n) — V(n)] s}. The coincidence holds 
also when importance sampling is included. If we change 
in Eqs. (|l6|-p"8"|) f2 into — ifl and interpret the absolute 
value as the modulus, in the continuum limit f2 -1 — > 
the extended GFMC algorithm coincides with our algo- 
rithm also for real times. 
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FIG. 1. Absolute value of the relative error 
(E — 75g xact ) /E^"" ct on the ground state energy of an inter- 
acting (7 = 4) and noninteracting (7 = 0) Hubbard system 
obtained with the present algorithm as a function of the Pois- 
son parameter p. The system is one-dimensional with peri- 
odic boundary conditions, |A| = 8, "■ x ^ 
rjij = 1 X Sij-i, and 74 = 4 or 74 = 0. 
used K = 10 3 , R = 10 4 and t = 0.26. 
statistical errors. 
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